clc;clear all;close all;
global alpha betta depr sigmaC sigmaL phi govexp omega taumax N T ALPHA popweight
global kstst lstst cstst rstst wstst kstart kistart tauK0
global Delta1 Delta2 lamda kg lam kmax

%run mainbench.m;

load('benchN2.mat')

kistart = kj;
kstart = sum(popweight.*kistart)+kg;
kmax=3;


lamda=lam;
Delta1 = .3;
Delta2 = -0.1;

kstst = fzero('findstst',kstart+0.3,options)

rstst = betta^(-1)-1+depr; %r from Euler at the steady state
estst = (rstst/alpha)^(1/(1-alpha))*kstst; % express total effective labor from r=F_k
for jj=2:N
ljpart(jj-1)=popweight(jj)*phi(jj)*lam(jj-1)^(-sigmaC/sigmaL)*(phi(jj)/phi(1))^(1/sigmaL);
end
lstst=estst/(popweight(1)*phi(1)+sum(ljpart)); %express l1 using the expression for total effective labor and l2=f(lam,l1)
for jj=2:N
    L2a(jj-1) = lamda(jj-1)^(-sigmaC/sigmaL)*(phi(jj)/phi(1))^(1/sigmaL); %last terms in (11) and f_l
end

T=150; %we assume that the economy reaches the steady state after T periods

tauK0 = taumax;

%initial guesses 
for i=1:T;
    time1 = i/T;
    X(i) = time1*(kstst) + (1-time1)*kstart; %k
    X(T+i) = lstst; %l
    if (N==5)
    X(2*T+i) = max(2-.07*i,0); %gamma
    end
    if (N==2)
    X(2*T+i) = max(1-.05*i,0); %gamma
    end
end;
X(3*T+1:3*T+1+2*(N-1)) = [Delta1 Delta2 lamda];

ALPHall=[1 0.95 0.9 0.8 0.7 0.6 0.5 0.488 0.487 0.4865 0.4863 0.4862 0.4861 0.48 0.46 0.43 0.4 0.37 0.35 0.347 0.3467 0.3466 0.3 0.25]
ii=1;    
%capitalist max U at 0.3467 (ii=21)
%worker max U at 0.4861 (ii=13)

for ii=1:length(ALPHall)

ALPHA=ALPHall(ii)
    
    if (ii>1)
    X=XX(ii-1,:);
    if ii==17
         X=0.5*XX(ii-1,:)+0.5*XX(ii-2,:);
    end
    Delta1 = X(3*T+1);
    Delta2 = X(3*T+2:3*T+2+N-2);
    lamda = X(3*T+3+N-2:3*T+3+2*(N-2));
    end
    for i=1:T;
        time1 = i/T;
        X(i) = time1*(kstst) + (1-time1)*kstart; %k
    end
kstst = fzero('findstst',kstart+0.3)
rstst = betta^(-1)-1+depr; %r from Euler at the steady state
estst = (rstst/alpha)^(1/(1-alpha))*kstst; % express total effective labor from r=F_k
for jj=2:N
ljpart(jj-1)=popweight(jj)*phi(jj)*lam(jj-1)^(-sigmaC/sigmaL)*(phi(jj)/phi(1))^(1/sigmaL);
end
lstst=estst/(popweight(1)*phi(1)+sum(ljpart)); %express l1 using the expression for total effective labor and l2=f(lam,l1)
for jj=2:N
    L2a(jj-1) = lamda(jj-1)^(-sigmaC/sigmaL)*(phi(jj)/phi(1))^(1/sigmaL); %last terms in (11) and f_l
end

tauK0 = taumax;
tolf0 = 1e-7;

if ii>=21

options = optimset('MaxFunEvals',10000000,'MaxIter',100,'TolFun',1e-15,'TolX',1e-15,'Disp','iter');
[X,check] = fsolve('resid3v',X,options)
X2=X;
X=X2;
if ii~=24
k = X(1:T);
l = X(T+1:2*T);
gamma = X(2*T+1:3*T);
Delta1 = X(3*T+1);
Delta2 = X(3*T+2:3*T+2+N-2);
lamda = X(3*T+3+N-2:3*T+3+2*(N-2));
kstst = fzero('findstst',[kstart-0.2 kmax],options)

for i=1:T;
    time1 = i/T;
    X(i) = time1*(kstst) + (1-time1)*kstart; %k
end
for i=1:T;
    X(2*T+i) = max(8-0.2*i,0); %gamma
end;

options = optimset('MaxFunEvals',10000000,'MaxIter',100,'TolFun',1e-15,'TolX',1e-15,'Disp','iter');
[X,check] = fsolve('resid3v',X,options)
end
end

if ii==20

options = optimset('MaxFunEvals',10000000,'MaxIter',100,'TolFun',1e-15,'TolX',1e-15,'Disp','iter');
[X,check] = fsolve('resid3v',X,options)

k = X(1:T);
l = X(T+1:2*T);
gamma = X(2*T+1:3*T);
Delta1 = X(3*T+1);
Delta2 = X(3*T+2:3*T+2+N-2);
lamda = X(3*T+3+N-2:3*T+3+2*(N-2));
kstst = fzero('findstst',[kstart-0.2 kmax],options)

for i=1:T;
    time1 = i/T;
    X(i) = time1*(kstst) + (1-time1)*kstart; %k
end
for i=1:T;
    X(2*T+i) = max(8-0.2*i,0); %gamma
end;

options = optimset('MaxFunEvals',10000000,'MaxIter',100,'TolFun',1e-15,'TolX',1e-15,'Disp','iter');
[X,check] = fsolve('resid3v',X,options)

end

if ii==19

for i=1:T;
    X(2*T+i) = max(5-0.1*i,0); %gamma
end;
options = optimset('MaxFunEvals',10000000,'MaxIter',100,'TolFun',1e-15,'TolX',1e-15,'Disp','iter');
[X,check] = fsolve('resid3v',X,options)

end


if ii==18
    
for i=1:T;
    X(2*T+i) = max(5-0.1*i,0); %gamma
end;
options = optimset('MaxFunEvals',10000000,'MaxIter',100,'TolFun',1e-15,'TolX',1e-15,'Disp','iter');
[X,check] = fsolve('resid3v',X,options)

for i=1:T;
    time1 = i/T;
    X(i) = time1*(kstst) + (1-time1)*kstart; %k
end
[X,check] = fsolve('resid3v',X,options)

end

if ii==16

for i=1:T;
    time1 = i/T;
    X(i) = time1*(kstst) + (1-time1)*kstart; %k
end
options = optimset('MaxFunEvals',10000000,'MaxIter',100,'TolFun',1e-15,'TolX',1e-15,'Disp','iter');
[X,check] = fsolve('resid3v',X,options)

for i=1:T;
    time1 = i/T;
    X(i) = time1*(kstst) + (1-time1)*kstart; %k
end
[X,check] = fsolve('resid3v',X,options)

save('solN2a.mat')
clc;clear all;close all;
load('solN2a.mat')


[X,check] = broydn('resid3v',X,tolf0)
if size(X,1)>1
    X=X';
end
sum(resid3v(X).^2)
[X,check] = fsolve('resid3v',X,options)
end

if ii==15

[X,check] = broydn('resid3v',X,tolf0)
if size(X,1)>1
    X=X';
end
sum(resid3v(X).^2)

for i=1:T;
    time1 = i/T;
    X(i) = time1*(kstst) + (1-time1)*kstart; %k
end

options = optimset('MaxFunEvals',10000000,'MaxIter',100,'TolFun',1e-15,'TolX',1e-15,'Disp','iter');
[X,check] = fsolve('resid3v',X,options)

end


if ii==14

options = optimset('MaxFunEvals',10000000,'MaxIter',100,'TolFun',1e-15,'TolX',1e-15,'Disp','iter');
[X,check] = fsolve('resid3v',X,options)
[X,check] = broydn('resid3v',X,tolf0)
if size(X,1)>1
    X=X';
end
sum(resid3v(X).^2)

options = optimset('MaxFunEvals',10000000,'MaxIter',100,'TolFun',1e-15,'TolX',1e-15,'Disp','iter');
[X,check] = fsolve('resid3v',X,options)
X2=X;
X=X2;
[X,check] = broydn('resid3v',X,tolf0)
if size(X,1)>1
    X=X';
end
sum(resid3v(X).^2)
for i=1:T;
    time1 = i/T;
    X(i) = time1*(kstst) + (1-time1)*kstart; %k
end

options = optimset('MaxFunEvals',10000000,'MaxIter',100,'TolFun',1e-15,'TolX',1e-15,'Disp','iter');
[X,check] = fsolve('resid3v',X,options)

end


if ((ii>=9) && (ii<=13))

options = optimset('MaxFunEvals',10000000,'MaxIter',100,'TolFun',1e-15,'TolX',1e-15,'Disp','iter');
[X,check] = fsolve('resid3v',X,options)


[X,check] = broydn('resid3v',X,tolf0)
if size(X,1)>1
    X=X';
end

options = optimset('MaxFunEvals',10000000,'MaxIter',100,'TolFun',1e-15,'TolX',1e-15,'Disp','iter');
[X,check] = fsolve('resid3v',X,options)

end


if ii==8

[X,check] = broydn('resid3v',X,tolf0)
if size(X,1)>1
    X=X';
end

options = optimset('MaxFunEvals',10000000,'MaxIter',100,'TolFun',1e-15,'TolX',1e-15,'Disp','iter');
[X,check] = fsolve('resid3v',X,options)

for i=1:T;
    time1 = i/T;
    X(i) = time1*(kstst) + (1-time1)*kstart; %k
end
for i=1:T;
    X(2*T+i) = X(2*T+i)+0.3; %gamma
end;
[X,check] = fsolve('resid3v',X,options)


for i=1:T;
    X(2*T+i) = X(2*T+i)+0.1; %gamma
end;
for i=1:T;
    time1 = i/T;
    X(i) = time1*(kstst) + (1-time1)*kstart; %k
end

[X,check] = fsolve('resid3v',X,options)

[X,check] = broydn('resid3v',X,tolf0)
if size(X,1)>1
    X=X';
end
[X,check] = fsolve('resid3v',X,options)

end

if ii==7

for i=1:T;
    time1 = i/T;
    X(i) = time1*(kstst) + (1-time1)*kstart; %k
end
options = optimset('MaxFunEvals',10000000,'MaxIter',100,'TolFun',1e-15,'TolX',1e-15,'Disp','iter');
[X,check] = fsolve('resid3v',X,options)


for i=1:T;
    time1 = i/T;
    X(i) = time1*(kstst) + (1-time1)*kstart; %k
end
[X,check] = fsolve('resid3v',X,options)
end


if ii==6
    
for i=1:T;
    time1 = i/T;
    X(i) = time1*(kstst) + (1-time1)*kstart; %k
end
options = optimset('MaxFunEvals',10000000,'MaxIter',100,'TolFun',1e-15,'TolX',1e-15,'Disp','iter');
[X,check] = fsolve('resid3v',X,options)

[X,check] = broydn('resid3v',X,tolf0)
if size(X,1)>1
    X=X';
end
sum(resid3v(X).^2)

for i=1:T;
    X(2*T+i) = max(5-0.1*i,0); %gamma
end;

[X,check] = broydn('resid3v',X,tolf0)
if size(X,1)>1
    X=X';
end
sum(resid3v(X).^2)

for i=1:T;
    time1 = i/T;
    X(i) = time1*(kstst) + (1-time1)*kstart; %k
end
options = optimset('MaxFunEvals',10000000,'MaxIter',100,'TolFun',1e-15,'TolX',1e-15,'Disp','iter');
[X,check] = fsolve('resid3v',X,options)

for i=1:T;
    time1 = i/T;
    X(i) = time1*(kstst) + (1-time1)*kstart; %k
end
for i=1:T;
    X(2*T+i) = max(5-0.1*i,0)+0.2; %gamma
end;
options = optimset('MaxFunEvals',10000000,'MaxIter',100,'TolFun',1e-15,'TolX',1e-15,'Disp','iter');
[X,check] = fsolve('resid3v',X,options)


[X,check] = broydn('resid3v',X,tolf0)
if size(X,1)>1
    X=X';
end
sum(resid3v(X).^2)

options = optimset('MaxFunEvals',10000000,'MaxIter',100,'TolFun',1e-15,'TolX',1e-15,'Disp','iter');
[X,check] = fsolve('resid3v',X,options)

end

if ii==5

options = optimset('MaxFunEvals',10000000,'MaxIter',100,'TolFun',1e-15,'TolX',1e-15,'Disp','iter');
[X,check] = fsolve('resid3v',X,options)


[X,check] = broydn('resid3v',X,tolf0)
if size(X,1)>1
    X=X';
end
sum(resid3v(X).^2)

options = optimset('MaxFunEvals',10000000,'MaxIter',100,'TolFun',1e-15,'TolX',1e-15,'Disp','iter');
[X,check] = fsolve('resid3v',X,options)
end


if ii==4

options = optimset('MaxFunEvals',10000000,'MaxIter',100,'TolFun',1e-15,'TolX',1e-15,'Disp','iter');
[X,check] = fsolve('resid3v',X,options)

for i=1:T;
    time1 = i/T;
    X(i) = time1*(kstst) + (1-time1)*kstart; %k
end
for i=1:T;
    X(2*T+i) = X(2*T+i)+0.3; %gamma
end;

[X,check] = broydn('resid3v',X,tolf0)
if size(X,1)>1
    X=X';
end
options = optimset('MaxFunEvals',10000000,'MaxIter',100,'TolFun',1e-15,'TolX',1e-15,'Disp','iter');
[X,check] = fsolve('resid3v',X,options)

end

if ii==3
    
k = X(1:T);
l = X(T+1:2*T);
gamma = X(2*T+1:3*T);
Delta1 = X(3*T+1);
Delta2 = X(3*T+2:3*T+2+N-2);
lamda = X(3*T+3+N-2:3*T+3+2*(N-2));
kstst = fzero('findstst',[kstart-0.2 kmax],options)

for i=1:T;
    time1 = i/T;
    X(i) = time1*(kstst) + (1-time1)*kstart; %k
end

options = optimset('MaxFunEvals',10000000,'MaxIter',100,'TolFun',1e-15,'TolX',1e-15,'Disp','iter');
[X,check] = fsolve('resid3v',X,options)


[X,check] = broydn('resid3v',X,tolf0)%,0,1);
if size(X,1)>1
    X=X';
end

for i=1:T;
    X(2*T+i) = max(6-0.1*i,0); %gamma
end;
for i=1:T;
    time1 = i/T;
    X(i) = time1*(kstst) + (1-time1)*kstart; %k
end
options = optimset('MaxFunEvals',10000000,'MaxIter',100,'TolFun',1e-15,'TolX',1e-15,'Disp','iter');
[X,check] = fsolve('resid3v',X,options)


for i=1:T;
    X(2*T+i) = max(6-0.1*i,0); %gamma
end;
for i=1:T;
    time1 = i/T;
    X(i) = time1*(kstst) + (1-time1)*kstart; %k
end

[X,check] = broydn('resid3v',X,tolf0)%,0,1);
if size(X,1)>1
    X=X';
end
sum(resid3v(X).^2)

options = optimset('MaxFunEvals',10000000,'MaxIter',100,'TolFun',1e-15,'TolX',1e-15,'Disp','iter');
[X,check] = fsolve('resid3v',X,options)
end


if ii==2
    
options = optimset('MaxFunEvals',10000000,'MaxIter',200,'TolFun',1e-15,'TolX',1e-15,'Disp','iter');
[X,check] = fsolve('resid3v',X,options)

for i=1:T;
    time1 = i/T;
    X(i) = time1*(kstst) + (1-time1)*kstart; %k
end
for i=1:T;
    X(2*T+i) = X(2*T+i)+0.3; %gamma
end;

options = optimset('MaxFunEvals',10000000,'MaxIter',100,'TolFun',1e-10,'TolX',1e-10,'Disp','iter');
[X,check] = fsolve('resid3v',X,options)

for i=1:T;
    X(2*T+i) = max(4-0.1*i,0); %gamma
end;
for i=1:T;
    time1 = i/T;
    X(i) = time1*(kstst) + (1-time1)*kstart; %k
end
[X,check] = fsolve('resid3v',X,options)

end

if ii==1 | ii==17
    options = optimset('MaxFunEvals',10000000,'MaxIter',200,'TolFun',1e-15,'TolX',1e-15,'Disp','iter');
    [X,check] = fsolve('resid3v',X,options)
end


XX(ii,:)=X;

%for ii=1:length(ALPHall)
%X=XX(ii,:);
k = X(1:T);
l = X(T+1:2*T);
gamma = X(2*T+1:3*T);
Delta1 = X(3*T+1);
Delta2 = X(3*T+2:3*T+2+N-2);
lamda = X(3*T+3+N-2:3*T+3+2*(N-2));
Delta1*kistart(1)+sum(Delta2.*kistart(2:N))
kstst = fzero('findstst',[kstart-0.2 kmax],options)%,0,0)


% save('solN2.mat')
% clc;clear all;close all;
% load('solN2.mat')

rstst = betta^(-1)-1+depr; %r from Euler at the steady state
estst = (rstst/alpha)^(1/(1-alpha))*kstst; % express total effective labor from r=F_k
for jj=2:N
ljpart(jj-1)=popweight(jj)*phi(jj)*lam(jj-1)^(-sigmaC/sigmaL)*(phi(jj)/phi(1))^(1/sigmaL);
end
lstst=estst/(popweight(1)*phi(1)+sum(ljpart)); %express l1 using the expression for total effective labor and l2=f(lam,l1)
for jj=2:N
    L2a(jj-1) = lamda(jj-1)^(-sigmaC/sigmaL)*(phi(jj)/phi(1))^(1/sigmaL); %last terms in (11) and f_l
end
wstst = (1-alpha)*kstst^alpha*estst^(-alpha); %w=F_e
cstst = 1/(popweight(1)+sum(popweight(2:N).*lamda))*(kstst.^alpha*estst.^(1-alpha) - depr*kstst - govexp); %c from resource constraint
mustst = omega*lstst^sigmaL*(1+sum(ALPHA.*lamda.^(-sigmaC).*phi(2:N)'/phi(1).*L2a)+(Delta1+sum(Delta2.*phi(2:N)'/phi(1).*L2a))*(1+sigmaL))/(wstst*(popweight(1)*phi(1)+sum(ljpart)));

klast = [kstart X(1:T-1)]; %k_{t-1}
e=[l' (l'*(lamda.^(-sigmaC/sigmaL).*(phi(2:N)/phi(1))'.^(1/sigmaL)))]*(popweight'.*phi);
e=e';

r = alpha*klast.^(alpha-1).*e.^(1-alpha); %r=F_k
w = (1-alpha)*klast.^(alpha).*e.^(-alpha); %w=F_e
c = 1/(popweight(1)+sum(popweight(2:N).*lamda))*(klast.^alpha.*e.^(1-alpha)+ (1-depr)*klast - k - govexp); %resource constraint

cnext = [c(2:end) cstst];
rnext = [r(2:end) rstst];
lnext = [X(T+2:2*T) lstst];

%check bound on tauk
constrtest = c.^(-sigmaC)./cnext.^(-sigmaC) - betta*(1+ (rnext-depr)*(1-taumax));

cistst = [cstst lamda*cstst]; %consumption of groups in the long run
listst = [lstst, L2(lamda,lstst,phi(2:N))']; %hours worked of groups in the long run
for tt=1:T
ci(tt,:) = [c(tt) lamda*c(tt)]; %consumption of groups during the transition
li(tt,:) = [l(tt), L2(lamda,l(tt),phi(2:N))']; %hours worked of groups during the transition
end

%lifetime utilities
if (sigmaC==1)
    V = sum(betta.^[0:T-1]'*ones(1,N).*(log(ci)- omega*li.^(1+sigmaL)/(1+sigmaL)))...
    + betta^(T)/(1-betta)*(log(cistst) - omega*listst.^(1+sigmaL)/(1+sigmaL));
else
    V = sum(betta.^[0:T-1]'*ones(1,N).*((ci.^(1-sigmaC))/(1-sigmaC) - omega*li.^(1+sigmaL)/(1+sigmaL)))...
    + betta^(T)/(1-betta)*((cistst.^(1-sigmaC))/(1-sigmaC) - omega*listst.^(1+sigmaL)/(1+sigmaL));
end

if sum(V>=Vsq)==2
    PI(ii)=1
else
    PI(ii)=0
end

residall(ii)=sum(resid3v(X).^2)
Vall(ii,:)=V
Vsq

%tax rates
tauKt = 1 - ((c(1:end-1)./c(2:end)).^(-sigmaC)/betta-1)./(r(2:end)-depr); %Euler
tauKt = [tauK0 tauKt];
plot(tauKt)
PI(ii)
tauKt
lamda_all(ii)=lamda;
Delta1_all(ii)=Delta1;
Delta2_all(ii)=Delta2;
term_all(ii)=1+ALPHA*lamda^(1-sigmaC)+(Delta1+Delta2*lamda)*(1-sigmaC);
kappa=lamda^(-sigmaC/sigmaL)*(phi(2)/phi(1))^(1/sigmaL);
terml_all(ii)=1+ALPHA*kappa^(1+sigmaL)+(Delta1+phi(2)/phi(1)*kappa*Delta2)*(1+sigmaL);
cst_all(ii)=cstst;
gammast_all(ii)=gamma(T);
must_all(ii)=mustst;

%foc for labor
klam=lamda^(-sigmaC/sigmaL)*(phi(2)/phi(1))^(1/sigmaL);
lamderl(ii)=-lamda^(-sigmaC/sigmaL)*(phi(2)/phi(1))^((1+sigmaL)/sigmaL)*betta^T*(1+sigmaL)*lstst^(sigmaL)/(sum(betta.^[0:T-1].*c.^(-sigmaC).*c)+betta^T/(1-betta)*cstst^(-sigmaC)*cstst...
    -sigmaC/sigmaL*lamda^(-(sigmaC+sigmaL)/sigmaL)*(phi(2)/phi(1))^((1+sigmaL)/sigmaL)*(sum(betta.^[0:T-1].*(-omega*l.^(sigmaL).*l))+betta^T/(1-betta)*(-omega)*lstst^sigmaL*lstst));
%foc for consumption
lamderc(ii)=-lamda*betta^T*(1-sigmaC)*cstst^(-sigmaC)/(sum(betta.^[0:T-1].*c.^(-sigmaC).*c)+betta^T/(1-betta)*cstst^(-sigmaC)*cstst...
    -sigmaC/sigmaL*lamda^(-(sigmaC+sigmaL)/sigmaL)*(phi(2)/phi(1))^((1+sigmaL)/sigmaL)*(sum(betta.^[0:T-1].*(-omega*l.^(sigmaL).*l))+betta^T/(1-betta)*(-omega)*lstst^sigmaL*lstst));
%end

save('solN2.mat')

end

clc;clear all;close all;
load('solN2_fb.mat')

term1=(1-sigmaC)*((1-betta)*Vcap_fb(18:36)+omega*listst_sq(1)^(1+sigmaL)/(1+sigmaL));
epscap_fb=((term1+1).^(1/(1-sigmaC))/cistst_sq(1)-1)*100
term2=(1-sigmaC)*((1-betta)*Vwor_fb(18:36)+omega*listst_sq(2)^(1+sigmaL)/(1+sigmaL));
epswor_fb=((term2+1).^(1/(1-sigmaC))/cistst_sq(2)-1)*100

load('solN2.mat')

term1=(1-sigmaC)*((1-betta)*Vcap(1:24)+omega*listst_sq(1)^(1+sigmaL)/(1+sigmaL));
epscap=((term1+1).^(1/(1-sigmaC))/cistst_sq(1)-1)*100
term2=(1-sigmaC)*((1-betta)*Vwor(1:24)+omega*listst_sq(2)^(1+sigmaL)/(1+sigmaL));
epswor=((term2+1).^(1/(1-sigmaC))/cistst_sq(2)-1)*100
epscap=[epscap0 epscap']
epswor=[epswor0 epswor']

plot(epswor(14:21),epscap(14:21),'-k','LineWidth',2.5);
axis([-1 7.6 -2.4 2.6])
hold on
plot(epswor(1:13),epscap(1:13),'-.k','LineWidth',2);
hold on
plot(epswor_fb(8:14),epscap_fb(8:14),'--k','LineWidth',1.5);
hold on
plot(epswor_fb(1:8),epscap_fb(1:8),'-.k','LineWidth',1);
hold on
plot(epswor(2),epscap(2),'-k','Marker','.','MarkerSize',25,'LineWidth',1.5);
hold on
plot(epswor(21:24),epscap(21:24),'-.k','LineWidth',2);
hold on
plot(epswor_fb(14:19),epscap_fb(14:19),'-.k','LineWidth',1);
hold on
plot([-1 7.6],[0 0],'--k');
hold on
plot([0 0],[-2.4 2.6],'--k');
xlabel("workers' welfare increase (percent)",'FontSize',12)
ylabel("capitalists' welfare increase (percent)",'FontSize',12)
legend('POPI','PO, not PI','with optimal lump-sum redistributive transfer PI','with optimal lump-sum redistributive transfer, not PI')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc;clear all;close all;
load('solN2.mat')

year=1:1:T;

X=XX(13,:);

k = X(1:T);
l = X(T+1:2*T);
gamma = X(2*T+1:3*T);
Delta1 = X(3*T+1);
Delta2 = X(3*T+2:3*T+2+N-2);
lamda = X(3*T+3+N-2:3*T+3+2*(N-2));

kstst = fzero('findstst',[kstart-0.2 kmax],options);%,0,0)
rstst = betta^(-1)-1+depr; %r from Euler at the steady state
estst = (rstst/alpha)^(1/(1-alpha))*kstst; % express total effective labor from r=F_k
for jj=2:N
ljpart(jj-1)=popweight(jj)*phi(jj)*lam(jj-1)^(-sigmaC/sigmaL)*(phi(jj)/phi(1))^(1/sigmaL);
end
lstst=estst/(popweight(1)*phi(1)+sum(ljpart)); %express l1 using the expression for total effective labor and l2=f(lam,l1)
for jj=2:N
    L2a(jj-1) = lamda(jj-1)^(-sigmaC/sigmaL)*(phi(jj)/phi(1))^(1/sigmaL); %last terms in (11) and f_l
end
wstst = (1-alpha)*kstst^alpha*estst^(-alpha); %w=F_e
cstst = 1/(popweight(1)+sum(popweight(2:N).*lamda))*(kstst.^alpha*estst.^(1-alpha) - depr*kstst - govexp); %c from resource constraint

klast = [kstart X(1:T-1)]; %k_{t-1}
e=[l' (l'*(lamda.^(-sigmaC/sigmaL).*(phi(2:N)/phi(1))'.^(1/sigmaL)))]*(popweight'.*phi);
e=e';
r = alpha*klast.^(alpha-1).*e.^(1-alpha); %r=F_k
w = (1-alpha)*klast.^(alpha).*e.^(-alpha); %w=F_e
c = 1/(popweight(1)+sum(popweight(2:N).*lamda))*(klast.^alpha.*e.^(1-alpha)+ (1-depr)*klast - k - govexp); %resource constraint

cistst = [cstst lamda*cstst]; %consumption of groups in the long run
listst = [lstst, L2(lamda,lstst,phi(2:N))']; %hours worked of groups in the long run
for tt=1:T
    ci(tt,:) = [c(tt) lamda*c(tt)]; %consumption of groups during the transition
    li(tt,:) = [l(tt), L2(lamda,l(tt),phi(2:N))']; %hours worked of groups during the transition
end

k1=k;
e1=e;
ccap1=ci(:,1);
cwor1=ci(:,2);
tauKt = 1 - ((c(1:end-1)./c(2:end)).^(-sigmaC)/betta-1)./(r(2:end)-depr); %Euler
tauKt1 = [tauK0 tauKt];
tauLt1 = 1 - omega/phi(1)*l.^sigmaL./(c.^(-sigmaC).*w); %consumption-leisure FOC
%tax revenues at each time t
tauKrev = tauKt1.*(r-depr).*[kstart k(1:end-1)];
tauLrev = tauLt1.*w.*(li*(popweight'.*phi))';
totaltaxrev=tauKrev+tauLrev;
deficit1=govexp - totaltaxrev;

X=XX(17,:);

k = X(1:T);
l = X(T+1:2*T);
gamma = X(2*T+1:3*T);
Delta1 = X(3*T+1);
Delta2 = X(3*T+2:3*T+2+N-2);
lamda = X(3*T+3+N-2:3*T+3+2*(N-2));

kstst = fzero('findstst',[kstart-0.2 kmax],options);%,0,0)
rstst = betta^(-1)-1+depr; %r from Euler at the steady state
estst = (rstst/alpha)^(1/(1-alpha))*kstst; % express total effective labor from r=F_k
for jj=2:N
ljpart(jj-1)=popweight(jj)*phi(jj)*lam(jj-1)^(-sigmaC/sigmaL)*(phi(jj)/phi(1))^(1/sigmaL);
end
lstst=estst/(popweight(1)*phi(1)+sum(ljpart)); %express l1 using the expression for total effective labor and l2=f(lam,l1)
for jj=2:N
    L2a(jj-1) = lamda(jj-1)^(-sigmaC/sigmaL)*(phi(jj)/phi(1))^(1/sigmaL); %last terms in (11) and f_l
end
wstst = (1-alpha)*kstst^alpha*estst^(-alpha); %w=F_e
cstst = 1/(popweight(1)+sum(popweight(2:N).*lamda))*(kstst.^alpha*estst.^(1-alpha) - depr*kstst - govexp); %c from resource constraint

klast = [kstart X(1:T-1)]; %k_{t-1}
e=[l' (l'*(lamda.^(-sigmaC/sigmaL).*(phi(2:N)/phi(1))'.^(1/sigmaL)))]*(popweight'.*phi);
e=e';
r = alpha*klast.^(alpha-1).*e.^(1-alpha); %r=F_k
w = (1-alpha)*klast.^(alpha).*e.^(-alpha); %w=F_e
c = 1/(popweight(1)+sum(popweight(2:N).*lamda))*(klast.^alpha.*e.^(1-alpha)+ (1-depr)*klast - k - govexp); %resource constraint

cistst = [cstst lamda*cstst]; %consumption of groups in the long run
listst = [lstst, L2(lamda,lstst,phi(2:N))']; %hours worked of groups in the long run
for tt=1:T
ci(tt,:) = [c(tt) lamda*c(tt)]; %consumption of groups during the transition
li(tt,:) = [l(tt), L2(lamda,l(tt),phi(2:N))']; %hours worked of groups during the transition
end

k2=k
e2=e;
ccap2=ci(:,1);
cwor2=ci(:,2);
tauKt = 1 - ((c(1:end-1)./c(2:end)).^(-sigmaC)/betta-1)./(r(2:end)-depr); %Euler
tauKt2 = [tauK0 tauKt];
tauLt2 = 1 - omega/phi(1)*l.^sigmaL./(c.^(-sigmaC).*w); %consumption-leisure FOC
%tax revenues at each time t
tauKrev = tauKt2.*(r-depr).*[kstart k(1:end-1)];
tauLrev = tauLt2.*w.*(li*(popweight'.*phi))';
totaltaxrev=tauKrev+tauLrev;
deficit2=govexp - totaltaxrev;


X=XX(21,:);

k = X(1:T);
l = X(T+1:2*T);
gamma = X(2*T+1:3*T);
Delta1 = X(3*T+1);
Delta2 = X(3*T+2:3*T+2+N-2);
lamda = X(3*T+3+N-2:3*T+3+2*(N-2));

kstst = fzero('findstst',[kstart-0.2 kmax],options);%,0,0)
rstst = betta^(-1)-1+depr; %r from Euler at the steady state
estst = (rstst/alpha)^(1/(1-alpha))*kstst; % express total effective labor from r=F_k
for jj=2:N
ljpart(jj-1)=popweight(jj)*phi(jj)*lam(jj-1)^(-sigmaC/sigmaL)*(phi(jj)/phi(1))^(1/sigmaL);
end
lstst=estst/(popweight(1)*phi(1)+sum(ljpart)); %express l1 using the expression for total effective labor and l2=f(lam,l1)
for jj=2:N
    L2a(jj-1) = lamda(jj-1)^(-sigmaC/sigmaL)*(phi(jj)/phi(1))^(1/sigmaL); %last terms in (11) and f_l
end
wstst = (1-alpha)*kstst^alpha*estst^(-alpha); %w=F_e
cstst = 1/(popweight(1)+sum(popweight(2:N).*lamda))*(kstst.^alpha*estst.^(1-alpha) - depr*kstst - govexp); %c from resource constraint

klast = [kstart X(1:T-1)]; %k_{t-1}
e=[l' (l'*(lamda.^(-sigmaC/sigmaL).*(phi(2:N)/phi(1))'.^(1/sigmaL)))]*(popweight'.*phi);
e=e';
r = alpha*klast.^(alpha-1).*e.^(1-alpha); %r=F_k
w = (1-alpha)*klast.^(alpha).*e.^(-alpha); %w=F_e
c = 1/(popweight(1)+sum(popweight(2:N).*lamda))*(klast.^alpha.*e.^(1-alpha)+ (1-depr)*klast - k - govexp); %resource constraint

cistst = [cstst lamda*cstst]; %consumption of groups in the long run
listst = [lstst, L2(lamda,lstst,phi(2:N))']; %hours worked of groups in the long run
for tt=1:T
ci(tt,:) = [c(tt) lamda*c(tt)]; %consumption of groups during the transition
li(tt,:) = [l(tt), L2(lamda,l(tt),phi(2:N))']; %hours worked of groups during the transition
end

k3=k
e3=e;
ccap3=ci(:,1);
cwor3=ci(:,2);
tauKt = 1 - ((c(1:end-1)./c(2:end)).^(-sigmaC)/betta-1)./(r(2:end)-depr); %Euler
tauKt3 = [tauK0 tauKt];
tauLt3 = 1 - omega/phi(1)*l.^sigmaL./(c.^(-sigmaC).*w); %consumption-leisure FOC
%tax revenues at each time t
tauKrev = tauKt3.*(r-depr).*[kstart k(1:end-1)];
tauLrev = tauLt3.*w.*(li*(popweight'.*phi))';
totaltaxrev=tauKrev+tauLrev;
deficit3=govexp - totaltaxrev;

figure
subplot(421)
plot(year,k1,'-k','LineWidth',1.5);
axis([0 150 1.2 1.9])
hold on
plot(year,k2,'--k','LineWidth',1.5);
hold on
plot(year,k3,'-.k','LineWidth',1.5);
%xlabel('time (years)')
title('Capital')
subplot(422)
plot(year,e1,'-k','LineWidth',1.5);
axis([0 150 0.22 0.28])
hold on
plot(year,e2,'--k','LineWidth',1.5);
hold on
plot(year,e3,'-.k','LineWidth',1.5);
%xlabel('time (years)')
title('Aggregate labor')
subplot(423)
plot(year,ccap1,'-k','LineWidth',1.5);
axis([0 150 0.33 0.45])
hold on
plot(year,ccap2,'--k','LineWidth',1.5);
hold on
plot(year,ccap3,'-.k','LineWidth',1.5);
%xlabel('time (years)')
title('Consumption of capitalists')
subplot(424)
plot(year,cwor1,'-k','LineWidth',1.5);
axis([0 150 0.17 0.22])
hold on
plot(year,cwor2,'--k','LineWidth',1.5);
hold on
plot(year,cwor3,'-.k','LineWidth',1.5);
%xlabel('time (years)')
title('Consumption of workers')
subplot(425)
plot(year,tauKt1,'-k','LineWidth',1.5);
axis([0 150 -0.1 0.5])
hold on
plot(year,tauKt2,'--k','LineWidth',1.5);
hold on
plot(year,tauKt3,'-.k','LineWidth',1.5);
%xlabel('time (years)')
title('Capital income tax')
subplot(426)
plot(year,tauLt1,'-k','LineWidth',1.5);
axis([0 150 -0.2 0.5])
hold on
plot(year,tauLt2,'--k','LineWidth',1.5);
hold on
plot(year,tauLt3,'-.k','LineWidth',1.5);
xlabel('time (years)')
title('Labour income tax')
subplot(427)
plot(year,deficit1,'-k','LineWidth',1.5);
axis([0 150 -0.1 0.12])
hold on
plot(year,deficit2,'--k','LineWidth',1.5);
hold on
plot(year,deficit3,'-.k','LineWidth',1.5);
xlabel('time (years)')
title('Government deficit')
legend('only workers gain', 'both gain', 'only capitalists gain')
